TUTORIAL 6 : Exercise 2

Author:Andrey V. Brukhno (andrey.brukhno{at}stfc.ac.uk)

Calculating PMF for two charged colloids in ionic cloud


The objective of this exercise is to use EE or WL methods for calculating PMF between two highly charged nano-particles surrounded by their counterions (in implicit solvent).

The simulation protocol to carry out PMF calculation along the separation between the centers of nano-particles is analogous to that described in the previous exercise (TUTORIAL 6 : Exercise 1).

Exercise 2.1: PMF in the case of monovalent ions (WL/RE scheme)

Change to directory tutorial6-2/FED_HSR10_Q20_EDL1, where you will find the input files for this exercise.

Here we represent two nano-particles (moltype ‘colloid’) and their counterions (moltype ‘ions’) as hard spheres (VDW type ‘hs’) bearing charges, which is reflected in the FIELD file:

Two charged nanoparticles and their counterions (hard spheres)
atoms 2
SRCP core  1.0  20.00
SRIN core  1.0  -1.00
vdw  3
SRCP core  SRCP core  hs  20.0  0.0000 0.000
SRCP core  SRIN core  hs  11.5  0.0000 0.000
SRIN core  SRIN core  hs  3.0   0.0000 0.000

Note that we have to use a rather large box and CUTOFF value in order to accommodate two spheres (‘SRCP’) of diameter 20 Angstrom (2 nm) plus some tens of small ions (‘SRIN’) of diameter 3 Angstrom.

As before, in CONTROL we specify the usage of parallel tempering replica-exchange (with a T-step of 50 K) combined with Wang-Landau biasing scheme along ‘com2’ order parameter:

use repexch 4  50.0  500

# FED block is started by 'use fed <flavor> <fed_freq>' directive
use fed generic 1
fed method WL   0.005  0.50  500000  3 #3 398 0.5 #0.7071068  100000

# FED order 'COM2' subsection for POMF acting between two centers of mass (COM)
fed order param com2     200  20.  40.  1
com  sampling correction 1
com1 mol 1 atom 1  # use COM1 for a set of molecules or atoms within molecules
com2 mol 1 atom 2  # use COM2 for a set of molecules or atoms within molecules
fed order done

# close the 'use' block

temperature          290.0    #K
epsilon              78.7
ewald precision      1.e-5

steps                8000000  #16000000
equilibration        0        #10000

Note how the implicit solvent model is specified by using the dielectric constant of water (78.7). In this case we perform simulation in an NVT (not NPT!) ensemble, so no directive for volume moves is included.

We have to carry out a relatively long simulation (8 million steps or more) to allow for satisfactory sampling over the entire range of COM separations. Later on we will see how this can be achieved more efficiently by running simulations in sub-ranges, also known as “windows”.

Do not forget that a parallel (replica exchange) job must be run via the queue system:

[tutorial6-2/FED_HSR10_Q20_EDL1]$ sbatch parallel.sub

After the job is successfully finished, using gnuplot you should obtain the following graph:

[tutorial6-2/FED_HSR10_Q20_EDL1]$ gnuplot
gnuplot> plot [x=20:40] [y=-1:25] 'FEDDAT.000_005' u 1:2 w l lc 'black' t "Replica 1, WL iter 5", \
'FEDDAT.001_005' u 1:2 t "Replica 2, WL iter 5" w l lc 'red', 'FEDDAT.002_005' u 1:2 w l lc 'green' t "Replica 3,\
 WL iter 5", 'FEDDAT.003_005' u 1:2 w l lc 'blue' t "Replica 4, WL iter 5"

As before, you should aslo check the corresponding probability distributions:

gnuplot> plot [x=20:40] [y=5000:15000] 'FEDDAT.000_005' u 1:3 w l lc 'black' t "Replica 1, WL iter 5", \
'FEDDAT.001_005' u 1:3 t "Replica 2, WL iter 5" w l lc 'red', 'FEDDAT.002_005' u 1:3 w l lc 'green' t "Replica 3,\
 WL iter 5", 'FEDDAT.003_005' u 1:3 w l lc 'blue' t "Replica 4, WL iter 5"

Is a refinement run necessary?

We have to options:

  • either increase the simulation length (steps) to 16 million and re-run the entire simulation once again,
  • or store away the already obtained results in a newly created subdirectory and use the latest FED data (FEDDAT.00?_005 -> FEDDAT.00?) and configurations (REVCON.00? -> CONFIG.00?) as input for another simulation of the same length as before.

In the latter case, do the following:

[tutorial6-2/FED_HSR10_Q20_EDL1]$ mkdir first-run/
[tutorial6-2/FED_HSR10_Q20_EDL1]$ mv H* FED* TR* R* S* P* O* first-run/
[tutorial6-2/FED_HSR10_Q20_EDL1]$ for i in {0..3}; do cp first-run/FEDDAT.00${i}_005 ./FEDDAT.00${i}; done
[tutorial6-2/FED_HSR10_Q20_EDL1]$ for i in {0..3}; do cp first-run/REVCON.00${i} ./CONFIG.00${i}; done

Remember: you have to instruct DL_MONTE to read FEDDAT.00? files as input by flipping the <power> in the FED specification from ‘1’ to ‘-1’ in CONTROL file:

fed order param com2     200  20.  40.  -1

Exercise 2.2: PMF in the case of divalent ions (WL/RE scheme)

Change to directory tutorial6-2/FED_HSR10_Q20_EDL2, where the input files are found.

Check the input files and perform FED-MC simulations.

Repeat the above protocol and see what you get!


Questions to answer:

  • Can you explain in physical terms the drastic change in the behaviour of nanoparticle net interaction - from pure repulsion to subtle attraction - upon
    1. replacing monovalent with divalent counterions?
    2. lowering temperature in the case with divalent counterions?
  • How increasing the charge on nanoparticles from 20 to 30 would affect the PMF in either case?

Next exercise

TUTORIAL 6 : Exercise 3 - Umbrella sampling in subranges (windows) with the use of WHAM method